Algorithm for a microfluidic assembly line 
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Microfluidic technology has revolutionized the control of flows at small scales giving rise to new 
possibilities for assembling complex structures on the microscale. We analyze different possible al- 
gorithms for assembling arbitrary structures, and demonstrate that a sequential assembly algorithm 
can manufacture arbitrary 3D structures from identical constituents. We illustrate the algorithm 
by showing that a modified Hele-Shaw cell with 7 controlled flowrates can be designed to construct 
the entire English alphabet from particles that irreversibly stick to each other. 



Developing novel methods for assembling complex 
structures from small particles has been the focus of 
much recent investigation [l] |2] . Traditional approaches 
have mostly revolved around designing selective interac- 
tions between constituent elements of the assembly. The 
structure is then assembled in the absence of any external 
control if thermal fluctuations can drive the system to its 
energetic ground state [TJ |3] , although non-trivial energy 
landscapes often render this approach challenging. 

Here we consider a different possibility for assembly on 
small scales, in which microfluidic flow control is used to 
steer and assemble small particles into structures of high 
complexity. The basic idea follows from the observation 
that if we could construct an arbitrary time dependent 
flow field v{x^t)^ then particles in the flow could be ad- 
vected along arbitrary paths and moved to arbitrary loca- 
tions at a fixed time. This apparently allows to construct 
any complex structure, with the individual components 
binding irreversibly upon contact. 

Of course, the flow field iT(x, t) cannot be arbitrary; it 
must conserve mass and momentum [4], 

V-iT=0, -Vp + /iV^i; + p6 = 0. (1) 

where p is the pressure, ja the liquid viscosity and 6 is a 
volumetric force. The question of what structures can be 
built thus hinges on what flow fields can be produced us- 
ing current technology [5 , and what are the limits for the 
possible structures that can form with such flow fields. 
Volumetric forces can be produced using e.g. magnetic 
fields or optical tweezers [7 ; alternatively, the flow 
can be generated by inlets specifying v at the boundary 
of a cell. 

We analyze the case where pressure inlets around the 
boundary force the flow (Fig. [l]), and 6 = 0. Fluid 
mechanical constraints prohibit simultaneous control of 
many particles with this device; however, we demonstrate 
that a sequential assembly algorithm allows the assem- 
bly of arbitrary structures. The device can be designed 
to manufacture the entire English alphabet using 7 con- 
trolled flow rates in two dimensions. 

Linear Response to Boundary Forcing. Consider N 
particles suspended in the flow domain Q with their in- 
stantaneous positions being Xj, j = 1, 2, . . . , A^. Let the 
flow be forced on the boundary of the flow cell by a pre- 
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FIG. 1: (Color online) (a) Schematic: Five particles are ad- 
vected by the flow field in a circular cell. The flow field, 
visualized by its streamlines, is set up by eleven flowrates im- 
posed on the boundary (arrows indicating the strength and 
direction), (b/c) Two trajectories transporting particles from 
a fixed initial to their desired final position with the required 
flowrates. The linear trajectory (b) requires flowrates almost 
2 orders of magnitude higher than the trajectory in (c). This 
optimized trajectory minimizes dissipation. Flowrates are 
given in units of the rate required to transport a single particle 
in one time unit across the cell. 



scribed velocity iT(x, t), x G dQ. The linearity of ([T]) 
implies the velocity of the suspended particles is linear 
in the boundary forcing [4 , i.e. 

^ • f 



2 



where Fj is the force acting on the j^^ particle, possibly 
due to non- hydro dynamic interparticle interactions. The 
response coefficients Kj and Rjk depend on the geometry 
of the flow cell [8 and can be computed numerically. If, 
the boundary forcing occurs through M discrete inlets 
of area AS*, located at with prescribed velocities Vk^ 
k = 1,2,..., M, then Eqn. Q can be written in the 
discrete form 

R{x) X = Mix) • / + F, (3) 

with X = [xi, X2, . . . , xat], flowrates / = 
[vi,V2, . . .,vm]^S and M_^.^ = Kj{xi,X2, • . .:r]V;6e)- 

Eqn. ^ is an instantaneous linear relation between 
the imposed flowrates and the particle velocities. This 
implies that prescribed particle trajectories may be 
realized by imposing suitable / obtained by inverting 
The feasibility of this method then hinges crucially upon 
invert ibility of M, which in general needs not be square. 

Condition for Controlling Particle Trajectories. A 
necessary condition for inverting M is that the number 
of independent controls exceed the number of degrees of 
freedom. With N particles in 3 dimensions, there must 
be at least M = 3N + 1 flow inlets; 3A^ inlets control 
particle degrees of freedom, with the additional inlet en- 
forcing volume conservation. Similarly in two dimensions 
at least 2N + 1 flow inlets are required. 

In general, these algebraic conditions are not sufficient 
for the practicality of this assembly method. We will see 
below that in practice, the flowrates required to indepen- 
dently steer large number of particles can be too large to 
be practical, owing to the poor conditioning of the matrix 

M. 

Hele-Shaw Flow as a Specific Example. To explore 
this in more detail, we consider a speciflc example. We 
consider the fluid mechanics of typical microfluidic de- 
vices [9] where the vertical gap thickness H is much 
smaller than the lateral extent. In Hele-Shaw approxima- 
tion [8] , the velocity proflle across the gap is a parabola, 
with the gap-averaged velocity u being proportional to 
the gradient of the pressure fleld p satisfying Laplace's 
equation. A particle at position Xk responds to the flow 
by moving at a speed proportional to the local fluid ve- 
locity, Xk = Pu{xk), with (3 depending on the particle 
size and shape [8j . This linear dependence of the particle 
velocity on the gap-averaged fluid velocity is a general 
consequence of the linearity of Stokes flow. The quanti- 
tative value for /3 can be calculated in two different limits: 
If the particle is much smaller than the gap, it is advected 
by the local velocity, so P depends on the location of the 
particle relative to the walls. Alternatively, if the par- 
ticle approximately spans the gap width, the separation 
of scale underlying the Hele-Shaw approximation breaks 
down near the particle. The factor P can now be calcu- 
lated by solving the Stokes equations close to the particle, 
matching the solution with the parabolic Hele-Shaw flow 



in the far fleld. 

We consider the device depicted in Fig. 1(a), a circular 
domain of radius a with flowrates fk prescribed at the 
boundary at M discrete inlets with positions Rk. The 
velocity fleld at any position x is then given by 

Thus, the matrix M in Eqn. (|3| may be constructed 
by combining the N position dependent matrices B_{xi) 
corresponding to each particle. 

We flrst test whether this flow cell will allow simulta- 
neous control of all N particles by boundary inlets. The 
flow rates scale inversely with the duration of the assem- 
bly so that the scale for the flux depends on the chosen 
timescale. We scale flowrate by the flux required to move 
a single particle across the cell F = iraHS/r^ where S is 
the width of the inlet and the duration r of the pro- 
cess. Fig. [TJb) shows that if we require the particles to 
be transported in straight lines at constant velocity from 
their initial position to their flnal position, the required 
flowrates reach up to 30 times this value. This happens 
because when two or more particles are moving towards 
each other a distance e apart, a strain rate fleld of e/e is 
required. If the approach velocity is constant, the strain 
diverges as e ^ 0, implying diverging flowrates. Thus, 
whenever particles are brought together at constant ve- 
locity, large fluxes are required. 

Optimized Trajectories. We can try to reduce the flow 
rates by choosing the particle trajectories and speeds con- 
necting the initial and the flnal states to minimize this 
effect. Intuitively, we can decrease the speed of the par- 
ticles as they approach each other to minimize the re- 
quired flowrates. To flnd out if this is sufficient we com- 
pute the optimal trajectories connecting the initial to the 
desired particle conflguration by flnding the trajectories 
that minimize dissipation. The dissipation rate is given 

by 

For the discretized forcing w can be expressed as a 
quadratic form in the flowrates w — P ■ D- f with metric 
D being analytically given in terms of the inlet positions. 
We minimize the dissipation and thereby the flowrates 
under the constraints that the dynamics move the parti- 
cles according to their equation of motion Eqn. (|3| from 
their chosen initial state to their flnal state. This re- 
quires that we flnd the trajectories that minimize the 
Lagrangian 

C = j\t{f^-^-l-X^ -[i-K-f] -76^-/} (6) 

where A and 7 are Lagrange multipliers to enforce the 
constraints and = (1, . . . , 1). A enforces the equation 
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of motion, whereas 7 requires that the fluxes / satisfy 
volume conservation. 

To minimize H we consider a smah perturbation of the 
flowrates in the direction of steepest descent, i.e. / 
/ + A/ with 



A/ 



— j =-e{2£-/ + Mt.A-^. (7) 



To flx the unknown flelds 7 and A we require the vari- 
ation of C with respect to x, A and 7 to vanish for the 
updated flowrates and trajectory x ^ x + Ax. Eval- 
uating the stationarity condition at linear order in the 
increments allows to eliminate 7 from Eqn. ([7|) and se- 
lects a unique solution of the adjoint equations 
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which are numerically solved together with the evolution 
of X using MATLAB. [13^ 

The minimization starts from the linear trajectory in 
Fig.jljb) along which we choose velocity variations damp- 
ing large peak flowrates. Panel (c) shows the same 5 
particles following an optimized trajectory for which the 
flowrates are reduced by more than an order of magni- 
tude, to the desired range. The reduction is due to choos- 
ing e/e to approach a flnite number as e ^ 0. Thus, op- 
timizing trajectories can lead to simultaneous control of 
a modest number of particles using boundary flowrates. 

Nonetheless, this approach does not scale to larger 
number of particles; e.g, we could not flnd trajectories 
allowing simultaneous control of 13 particles required to 
spell the letter 'B' at moderate fluxes. This reflects a 
implementation independent physical limit of directing 
flow flelds with boundary fluxes, resulting from the fact 
that flow modes forced at the boundary decay in ampli- 
tude as one moves away from the wall. The character- 
istic length scale governing the decay of the boundary 
modes is the distance between the injection points which 
scales inversely with the number of inlets. But as the 
particle number increases, we need more inlets to control 
the particle motion. This limits the number of particles 
that can be controlled simultaneously. In our numerical 
experiments we could not flnd trajectories transporting 
more than 6 individual particles at moderate flowrates. 

Sequential Assembly. To overcome this physical limi- 
tation, an algorithm is required that decouples the num- 
ber of controlled degrees of freedom from the num- 
ber of particles in the desired structure. This can be 
achieved with a sequential approach, where one parti- 
cle after the other is attached to an aggregate which 
moves as a rigid body subject to force- and torque bal- 
ance. Construction of the corresponding M matrix re- 
quires explicitly accounting for the translation and rota- 
tion of the cluster consisting of Z particles at positions 

-R{^) {ia — icrn^ whcrc ^ are the prescribed 



particle positions in the aggregate's frame of reference 
and R is the two-dimensional rotation. How many de- 
grees of freedom are required for the assembly? Control- 
ling the position of the aggregate requires two degrees 
of freedom. Rotating it requires a stagnation point flow 
superimposed on the local mean flow near the aggregate. 
This requires independent control of two more degrees of 
freedom, the strength and orientation of the stagnation 
point flow. Finally, the free particle location demands 
two additional degrees of freedom. Including volume con- 
servation we thus need 7 inlets to absolutely control all 
degrees of freedom. 

Sequentially adding particles allows to assemble struc- 
tures of arbitrary shape and thus spell a word (Fig. [2]) 
without any feedback control. Technically, the trajectory 
for a particle to be added is constructed by deflning the 
desired position which is taken from [10 and the direction 
of approach in the frame of reference of the aggregate. 
Spline interpolation between the initial position of the 
new particle and its flnal position when the aggregate has 
reached the desired conflguration, then yields the particle 
path. The aggregate can also be moved arbitrarily; we 
choose to keep its orientation flxed and move the center 
of mass to the device center in each step. The velocity 
along this trajectory is chosen such that it vanishes at the 









FIG. 2: (Color online) Sequentially particles are added 
to form aggregates of arbitrary shape. The three columns 
present snapshots along the assembly of three example struc- 
tures presented in the last row. The required time-variation 
of the seven controlled flowrates are computed a priori (eg. 
Fig. [3]). See movie Ml (online) for an animation. 
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FIG. 3: (Color online) Flowrates fk/7vaH6T~'^ for the seven 
inlets as a function of time t/r. These allow to spell the letter 
'B'. At each integer time a new particles enters the cell. 



plement feedback control [12], though this significantly 
complicates the process especially for three dimensional 
structures. Another intriguing possibility is to simply 
embrace the existence of stochasticity and construct the 
most stable trajectories that maximize the probability 
of formation of the desired structures in the presence of 
noise. Finding algorithmic methods for carrying out this 
optimization are important directions for future research. 
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initial and final configuration so that flowrates reach zero 
after each assembly step. With this choice of trajectories, 
the flowrates are computed inverting Eqn. ([3|. No ex- 
tensive optimization was required to limit the fluxes (see 



3|) to reasonable values of the same order as those in 
T 



Fig. 

Fig. 1 'c). Optimizing trajectories might however still be 
useful either to further reduce the required flowrates or 
to minimize internal forces within the growing aggregate 
so that the forming structure can sustain the stresses in 
the flow. 

We envision a typical mode of operation of the mi- 
crofluidic device similar to that of a macroscopic robotic 
assembly line: Once the constituents of an assembly 
and the detailed assembly sequence are decided on, the 
flowrates on the inlets of the device and the precise in- 
stants at which each constituent is introduced are cal- 
culated. These pre-computed flowrates are then set up 
in a time-periodic fashion, and a train of the assembly 
constituents are introduced at the inlets of the device to 
accomplish a train of assembled products. 

In summary, we have shown that microfluidic assem- 
bly can be an efficient strategy if structures are built 
sequentially. Using this approach, we show that the tem- 
poral control of only 7 flowrates allows to build arbi- 
trarily shaped particle aggregates in 2 dimensions. In 
3 dimensions, the same argument implies that 11 dif- 
ferent flowrates are required. Different chamber geome- 
tries and hydrodynamic interactions can be incorporated 
into this framework, which rests solely on the linearity 
of Stokes flow. We anticipate that similar algorithms 
can be constructed using other forcing mechanisms such 
as electrokinetics [11 , where non- hydrodynamic electri- 
cal forcing contributions need to be included into the 
transfer matrices. A challenge for practical implementa- 
tion is to quantify the sensitivity of particle trajectories 
to noise in the imposed flow rates. This is particularly 
relevant when scaling down the assembly to submicron 
scale, where the Peclet number corresponding to Brow- 
nian motion is small and hence potentially disruptive. 
One option for dealing with errors and noise is to im- 
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